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Abstract 

Numerical  solutions  of  the  two-dimensional  averaged  Navier-Stokes 
equations  for  the  prediction  of  laminar,  transitional,  and  turbulent  flow 
fields  around  finite  bodies  of  arbitrary  shapes  have  been  considered. 
Numerically  generated  body-fitted  curvilinear  coordinates  and  the  relevant 
metric  terms  are  used  to  provide  the  finite-difference  solutions  of  the 
Navier-Stokes  and  the  equations  of  turbulent  quantities.  Complete  flow 
fields,  including  the  boundary  layer  parameters,  are  obtained  by  using 
the  zero,  one  and  two-equations  models  for  Schubauer's  elliptical  section, 
NACA66^-018  airfoil,  and  a  circular  cylinder  at  free  stream  Reynolds  numbers 
of  1.59  x  105,  1.2  x  106  and  1.4  x  106  respectively.  In  addition,  a  two- 
equation  model  with  an  algebraic-stress  closure  has  also  been  developed. 
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Function,  Eq.  (48) 

Constant 

Pressure  and  Skin  Friction  Coefficients 

Diffusion,  Eq.  (43) 

Differential  Operator,  Eq.  (lb) 

Turbulence  Energy  per  Unit  Mass 

e/u2 

00 

det  (g^) 

Function,  Eq.  (40d) 

Metric  Coefficients 
Function,  Eq.  (40c) 

0.41 

Characteristic  Length 
Mixing  Length 
Functions,  Eq.  (53) 

Normal  Distance  from  the  Surface 
n/L 

Pressure 

Production,  Eq.  (38) 

Function,  Eq.  (25) 

Function,  Eq.  (34) 

Function,  Eq.  (25) 

U  L 
00 

Reynolds  Number  — 

Turbulence  Reynolds  Number,  Eq.  (40) 
Reference  Value  of  R 
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-  R  7(1  +  R  T) 

o  o 

■  Function,  Eq.  (25) 

-  Length  Along  the  Surface 

■  Rate-of-Strain  Tensor 

*  Dimensional  and  Non-Dimensional  Times,  Respectively 

*  Eddy  Viscosity,  Eq.  (28) 

=  Vector,  Eq.  (18) 

■  Functions,  Eq.  (53) 

■  Velocity  Vector 

-  Non-Dimensional  Cartesian  Velocity  Components 

■  Free  Stream  Velocity 

■  Friction  Velocity 

■  Velocity  Along  the  Surface  Coordinates 

=  Max  (U  ) 
s 

*  Acceleration  Parameter 

=  Subscripted  Cartesian  Coordinates 

*  Non-Dimensional  Cartesian  Coordinates 

=  Non-Dimensional  Cartesian  Coordinates  in  Subscript  Form 
=  Roughness  Height 

*  Constants 

=  General  Coordinates 
=  General  Function 
=  Laplacian  in  x,y 
=  Laplacian  in  x^ 

-  Christoff el  Symbol,  (see  Appendix  A) 

■  Stream  Function 

■  Geometrical  Coefficients,  Eqs.  (9), (10) 


Tl*al  =  Modified  t  and  a  Respectively,  Defined  in  the  Text 

*  Dissipation,  Eq.  (25) 

X  *  Dissipation,  Eq.  (42) 

=  Diffusion,  Eq.  (25) 

=  Non-Dimensional  Reynolds  Stress 

a)  *  Vorticity 

“  -r 

w 

u)w  =  Wall  Vorticity 

v  =  Kinematic  Viscosity 

vT  =  Eddy  Viscosity 

6*, A*  -  Displacement  Thickness 

A  =  i/L 

Y*.  *  Transition  Factor 

tr 

Subscripts  and  Superscripts: 

(  .  )  *  Averaged  and  Dimensional  Quantity 

(  .  )"*  *  Perturbation 

(  7  )  =  Tensor 

(  .  )  =  Vector 

(  .  )^  =  Tensor  Components 
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1.  Introduction 

The  computation  of  mean  turbulent  flow  fields  for  flows  past  finite 
bodies  without  effecting  the  thin  shear  layer  approximation  is  a  problem  of 
much  computational  and  practical  importance.  There  are  three  distinct  areas 
of  immediate  attention  if  one  is  to  solve  the  problem  of  mean  turbulent 
flows.  The  first  is  the  problem  of  and  the  order  of  closure  which  can  be 
considered  as  satisfactory  for  a  given  problem.  The  second  is  in  regard  to 
the  generality  and  ease  with  which  the  body  geometry  can  be  brought  into 
the  process  of  calculations  so  that  arbitrary  shaped  bodies  can  be  considered 
as  simply  as  those  having  simpler  geometries.  The  third,  and  by  no  means 
lesser  than  the  above  two,  is  the  development  of  efficient  numerical 
techniques  which  can  give  accurate  solutions  in  reasonable  amounts  of  time. 

Most  of  our  present  state  of  knowledge  regarding  closure  is  usually  in 
the  context  of  boundary  layers,  far  wakes  and  jets.  Not  much  is  available 
by  the  way  of  experimental  data  for  turbulent  flows  in  the  separated  and 
recirculating  regions  near  a  body  surface.  Thus  the  modeling  of  various 
terms  which  have  proved  useful  in  the  context  of  boundary  layers  are 
heuristically  extended  to  apply  for  non-boundary  layer  regions  too.  How 
much  error  is  introduced  in  the  predicted  values  of  the  field  variables 
this  way  remains  to  be  checked  in  the  future.  A  recent  article  by  Reynolds 

[1]  summarizes  the  present  state-of-the-art  in  the  closure  problem. 

The  problem  of  arbitrary  shaped  bodies  and  general  coordinate  systems 
for  use  in  the  numerical  computations  of  the  viscous  flow  problems  is  not 
as  compelling  as  it  was  before  the  work  of  Thompson,  Thames,  and  Mastin 

[2] .  The  so  called  body-fitted  coordinates  method  developed  in  Ref.  [2]  is 
the  method  in  which  the  geometry  of  the  body  is  entirely  preserved  by  taking 
the  curve  forming  the  body  surface  as  one  of  the  coordinate  curves.  This 
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method  has  successfully  been  used  in  solving  various  fluid  dynamical 
problems  including  the  problem  of  supersonic  flow  past  a  blunt  body,  [3]. 
Though  at  present  various  authors  are  developing  other  techniques,  e.g., 
Shamroth,  et.  al.,  [4],  Eiseman  [5],  and  Davis  [6],  the  "differential 
equation"  approach  initiated  in  Ref.  [2]  looks  to  have  more  potential  for 
future  research  and  development. 

In  the  area  of  finite-difference  method  for  the  solution  of  fluid 
dynamical  equations  there  are  various  explicit  and  implicit  methods  of 
solution  available  in  the  open  literature,  e.g.,  Roache  [7].  The  choice 
of  a  method  usually  depends  on  the  formulation  of  the  overall  problem 
and  the  character  of  the  equations. 

The  main  aim  of  this  research  has  been  to  develop  a  suitable  numerical 
program  of  solution  for  the  solution  of  two-dimensional  mean  Navier-Stokes 
equations  in  non-orthogonal  curvilinear  coordinates.  The  coordinates  and 
all  the  metric  coefficients  are  assumed  to  be  available  in  numerical  form. 
The  body-fitted  coordinate  system  generation  technique  of  Ref.  [2]  has  been 
used  to  generate  coordinates  for  various  two-dimensional  shapes.  Care 
has  been  taken  to  have  the  maximum  concentration  of  coordinate  lines  in 
the  neighborhood  of  the  body  surface. 

In  the  present  research  three  body  shapes  have  been  considered  for 
the  prediction  of  laminar  transitional  and  turbulent  regions  using  different 
orders  of  modeling.  The  first-shape  is  an  elliptical  section  as  considered 
by  Schubauer  [8]  at  a  free  stream  Reynolds  number  of  158860.  Both  the  zero- 
and  one-equation  models  have  been  used  with  the  transitional  model  at  the 
predicted  points  on  the  body  surface. 

The  second  shape  is  the  NACA662-OI8  airfoil  at  zero  incidence  and  the 
chord  Reynolds  number  of  1.2  x  106.  In  this  case  also  the  zero-and  one- 
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equation  models  have  been  used.  The  energy  equation  has  been  modeled 
according  to  Glushko  [9]  in  both  these  problems. 

The  third  shape  is  a  circle  at  a  Reynolds  number  of  1.4  *  106  based 
on  the  radius.  In  this  problem  the  two-equations  model  of  Saffman  and 
Kolmogorov  [10]  has  been  used. 

Details  of  turbulence  modelings  are  given  in  Section  3  of  this  paper. 
In  the  course  of  our  investigation,  we  have  also  developed  a  method  of  two- 
equations  model  with  the  Reynolds  stress  algebraic  closure  which  can  be 
called  as  equivalent  to  solving  all  the  relevant  Reynolds  stress  equations. 
However,  because  of  the  storage  requirements  needed  by  this  method,  the 
method  could  not  be  tested  for  all  the  problems. 
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2.  Analysis 

This  section  briefly  describes  the  basic  aspects  of  numerical  coordinate 
generation  and  the  transformation  of  the  fluid  dynamical  equations  to  non- 
orthogonal  coordinate  systems.  For  a  description  of  th'j  quantities 
appearing  in  the  following  equations,  refer  to  the  list  of  symbols  and  to 
Appendix  A. 

All  the  pertinent  formulae  used  in  this  paper  stem  from  the  formula 
for  the  Laplacian  of  a  scalar  where 

<Kx,y)  =  <Kx(C,n),  y(£,n)). 


Here  x  and  y  are  the  Cartesian  coordinates  and  £x  =  £,  £2  =  n  are  the 
curvilinear  coordinates.  The  Laplacian  V2<(>  in  the  (£J,£2)  coordinates  becomes 


in  (la) ,  we  get 


V2<t»  =  i  D2(j>  -  r* 

g  9£r 


Setting  in  turn  $  =  £*  =  £  and  41  *  £2 


n  in  (lc),  we  get 


Again  setting  in  turn  $  *  x  and  $  ■  y  in  (Id) ,  we  get 


D2x  -  -g[(V2£)x5  +  (v2n)xnl  ) 

D2y  -  -g[(v2£)y  +  (v2n)y  ]  t  (2) 

£  n 

where  a  variable  subscript  denotes  a  partial  derivative. 

The  method  of  Thompson ,  Thames,  and  Mastin  [ 2  ]  for  the  two-dimensional 
numerical  coordinate  generation  now  depends  on  the  assumption  that  the 
coordinates  £  and  n  satisfy  an  elliptic  system  of  equations,  i.e., 

v2£  =  -fx(£,n)  ,  v2n  =  -f2(£,n)  ,  (3) 


where  the  functions  ^  and  f 2  are  treated  as  control  functions  to  be 
specified  in  a  manner  so  as  to  achieve  a  desired  distribution  of  the 
coordinates.  Thus  Eqs.  (2)  become 


D2x  =■ 

!  + 

f2xn) 

(4) 

D2y  = 

!  8^!^  + 

f2yn) 

(5) 

The  boundary  conditions  for  Eqs. 

(4)  and 

(5)  are 

on  T^*:  x  = 

Fi«-V 

,  y-F2«,nB)  | 

_  * 

on  Tj  :  x  = 

.  V  -  G2({,n_>  1 

(6) 

where  and  r£*  are  the  images  of  the  body  and  the  external  body  contours 
in  the  transformed  plane,  cf.  Figure  1. 


The  method  of  solution  of  Eqs.  (4)  and  (5)  under  the  boundary  conditions 
(6)  with  periodic  conditions  on  the  cut  has  been  described  in  various  reports, 
e.g.,  [llJ»  tl2 ] •  [13 ] >  and  will  not  be  described  here. 


Another  form  of  the  Laplacian  V2$ ,  to  be  used  in  the  transport  equations 


we  get 


where 


’  from  Eqs. 

(2).  Solving  Eqs. 

(2)  for  V2£  and  V2n 

=  x/g 

(7) 

v2n 

■  o/g 

(8) 

t  =  (x^D2y 

-  y^D2x)/g1^2 

(9) 

o  =  (y^D2x 

-  x^D2y) / g1/2 

(10) 

V2$  =  (d2<j>  +  +  o4>r))/g  (11) 

In  this  paper  we  have  used  the  vorticity-stream  function  formulation 
for  the  equations  of  motion,  which,  along  with  the  transport  equations  for 
turbulence  quantities  can  be  represented  as 

X  +ifiA  -  \p  X  =  V2X  +  B 
t  y  x  x  y 

where  all  quantifies  in  the  above  equation  are  non-dimensional;  p  is  the 
stream  function,  X  and  X  are  surrogate  variables  and  B  is  a  function  of 
the  dependent  variables.  Using  Eqs.  (7)  and  (8),  the  representative 
equation  transforms  as 

\  +  ~T/2  (Vs  ‘  Vn}  “  i(°2*  +  tXC  +  °V  +  B  ’  (12) 


; 


3.  Mean  Navier-Stokes  Equations  and  Turbulence  Modeling 

The  pertinent  equations  for  the  description  of  mean  turbulent  flow  are 
obtained  by  replacing  the  instantaneous  dependent  variables  in  the  continuity, 
Navier-Stokes  and  the  vorticity  equations  by  the  sum  of  a  mean  and  a 
perturbation  and  then  performing  the  time  average  of  each  term  of  the 
equations.  The  time  average  is  assumed  to  be  performed  to  remove  the 
turbulent  fluctuations  while  maintaining  the  time  dependency  of  the  mean 
flow.  Denoting  the  averages  by  an  overhead  bar  and  the  fluctuations  by 
a  prime,  we  introduce 

u  ■  u  +  u"  »p*p  +  p””,  u  ■  u>  +  w'',  etc.  (13) 

in  the  instantaneous  equations  and  after  averaging  obtain  the  following 
equations. 

Continuity: 

div  u  =  o  ,  div  ui  ®  o  (14) 

Momentum: 

3u  _  i  _  _  a 

-r-—  +  div  (uu)  - - grad  p  +  vV^u  -  div  t  (15a) 

ot  _  D  1" 


3u 

— —  +  div  (uu) 
3tl 


-  j  grad 


p  +  W^u  -  grad  e  +  T 


(15b) 


Vorticity: 

3u>  _  _  _  _  ___  _ 

-r-—  +  (grad  u>)  •  u  -  (grad  u)  *  to  ■  vVfu  +  div  (u'u'  -  w'u')  (16) 

-  *  -  l” 


where 


__li_ 

l 


e  is  the  mean  turbulence  energy  per  unit  mass. 
Comparing  (15a)  and  (15b) ,  we  get 

T  -  div  (el  -  t)  , 


(17) 


(18) 


where  I  is  the  unit  tensor.  Using  the  vector  identity 


we  have 


curl  (u'  x  u')  m  div  (u'cu'  -  oj'u') 


curl  T  ”  div  (u'w"  -  w'u')  , 


which  establishes  the  form  of  the  last  term  on  the  right  hand  side  of  (16). 

In  all  cases  studied  in  this  paper,  except  the  one  to  be  described  in 
3.1.4,  we  have  used  the  Newtonian  constitutive  equation  to  model  the 


Reynolds  stresses,  viz.. 


x  •  ^  el  -  2vtD 


where 


D  -  -|[grad  u  +  (grad  u)*] 


3.1  Transport  Equations  in  Cartesian  Coordinates 

In  the  case  of  Cartesian  coordinates,  x^  x^,  Xy  the  equations  are 


3u>.  _  3u> .  _  3u^  3T^ 

3t[  +  Uk  3^  '  “k  3^  "  v7l“i  +  eijk  3^ 


,  i  =  1,2,3  (21) 


^  +  \^;-5  +  6-I  +  v7ii 


+  u,  *  S  +  vV^e 

3^  k  3x^  1 


+  +  V  +  vV^13  1  11:1  '  1>2’3  '  W) 


where 


pij  "  -(,ik  a^  +  Tjk  to? 


8 


_  *•  3u? 

QU  ■  V  <8^  +  5^> 


'« ■  ■  si  lTn“k +  V  <u?ik + 


3u:  3U; 

e  =  2v  — -  — 1 
ij  9Xk  3Xk 


P  =  i  P 

2  ii 


D  -  Hi 


-  1  - 

E  "  2  eii 


3^u  3u;  3u  3u:  3uT 

-2v  [- — (u"  — +  — J-  ( — J - — 

3V*k  k  3X£  3\  <3x*  3Y 


3G.  3u;  3U;  3U;  3U;  9u; 

+  — t.  ( — 1 _ l-v  +  _i _ i _ k 

8xi  *xt  v  a*t  9-i  9xt 


i  «  3uf  o  ,  -  3u "  7 


+  V<8W  1 

Equations  (21)- (25)  form  a  closed  system  of  equations  when  the  terms 
occuring  in  (25)  have  been  modeled  on  the  basis  of  empirical— experimental 
and/or  heuristic  reasoning. 

In  this  paper,  we  have  considered  only  the  two-dimensional  mean  turbulent 
flows,  and,  except  for  one  case  where  Eq.  (24)  has  been  used  to  model  the 


Reynolds  stresses  ,  all  cases  use  the  concept  of  eddy  viscosity,  viz., 
Eqs.  (19),  (22),  and  (23).  For  two-dimensional  flows,  Eq.  (21)  is  written 
as 


^ + \  “  V  5  +  5 


where 


fl  -  (v  +  vT)Z 


B  -  2( 


.•N 

3u^ 

92vt 

_ £  4.  0 

•S 

% 

3x2 

3x^2 

3x2 

3x22 

-  T  Z. 

3xl 

3x^3x2 

where  v^,  is  the  kinematic  eddy  viscosity. 

We  now  introduce  the  following  scheme  for  non-dimensionalizing 
Eq.  (26): 

8  =  —  ,  U)  =  —  ,  T  =  — y  ,  <l>  =  -4r-  \ 

2  u  uL  uL  I 

U  oo  oo  oo  I 


U»C1  U1  u2 
c  =—  ’  u  mT  •  v  “IT 


X1  y!  p  U”L  /  (2R) 

x  '  T  ’  y  =  T  ’  Ro  '  “  /  (28) 

Where  u  is  the  freestream  mean  velocity  and  L  a  characteristic  length. 

00 

Thus  the  vorticity-stream  function  formulation  for  the  flow  field  in  terms 
of  the  non-dimensional  coordinates  takes  the  form  (cf.  Eq.  (12)) 

D2ip  +  +  ot|^  =  -goi  (29) 


“t  +  ~ffi  ‘  *(“»’  '  s(D2°  +  T“?  +  °“n>  +  * 
g 


where 


(l) 

amr 

U> 
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«  1  +  R  I 

o 


B  -  2(T  u  -  T  v  +  2T  v  ) 
xx  y  yy  x  xy  y 


and  the  variable  subscripts  S,n,x  and  y  denote  partial  derivatives. 
The  boundary  conditions  for  Eqs.  (29)  and  (30)  are: 


n  »  Hgt  m  0  *  'j'y  ■  0;  T  ■  0;  u)  to  be  calculated 
iteratively  to  satisfy  ^  ■  0  . 


n  -*•  n„:  i|»  -*■  y;  T-*- 0;  u»  -*•  0  , 


n  denote  the  wall  and  external  boundaries  respectively. 

B  * 

Below  we  now  consider  the  zero,  one  and  two-equation  models  and  also 
the  "algebraic-stress  closure"  model,  respectively. 


3.1.1  Zero-Equation  Modeling 

In  the  zero— equation  model,  only  Eqs.  (29)  and  (30)  are  solved  by 
specifying  the  eddy  viscosity  vT  by  a  generalized  mixing  length  model 
(cf.  [1],  [14])  and  defined  as 


(Vi  - *2(2  V 


(32a) 


(Vo 


.0168 | u  | 6* 


1  +  5.5(-S-)6 
n6 


(32b) 


where  l  is  the  mixing  length,  the  subscripts  i  and  o  denote  the  inner  and 
outer  values,  n  is  the  normal  distance  from  the  body,  ng  is  the  thickness 
of  the  shear  layer,  u^  is  the  velocity  tangential  to  the  surface  at  n  *  n^. 


1  3ui  3u1 
S  ■  — ( — -  +  — ■*•) 
ij  2l3Xj  3x^ 
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r 


r- 

F. 

F  ! 


n.  u 

6*  -  /  6(1  -  -®)dn 
o  - 

u 


u  *  max  (u  ) 
P  s 


Introducing  the  non-dimensional  quantities 


#  ■  r  - Q  ■  h <2  5«  V  - « - !  •  "s  •  t  ’ 


u  u  fit 

U  ■  —  ,  U  »  ,  A*  =  -y- 

8  “  P  u«  L 


We  have 


I  < 


In  two  dimensions 


(T)i  =■  A2Q1/2 

.0168 lu  I  A* 

<T)„ - Lt— 


M  6 

1  + 

<5 


Q  ■  (i<»  -i*<  )2  +  4  i|i 

H  xx  yy  xy 


oj2  +  4(i|/  2  -  ip  t  ) 
xy  xx  yy 


0  - 
8  F  g 


&12 

s  ♦,  -  JE=*< 


'88 


11 


N  U 

4*  *  L  (1  -  u£)<1" 

P 


(33) 


(34) 


(35) 


(35) 


The  mixing  length  is  given  by  the  Van  Driest  formula  for  the  near  wall 
region  Joined  with  the  Bradshaw's  [15]  empirical  fit  as  follows. 
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i _ I _ > 


A  -  kN[l  -  exp(-N^Ro|a)w|/26.)J  ,  o  <  N  <  .22Nj 

-  .089N6  tanh(4.6067N/N6)  ,  .22Nfi  <  N  <  •  (36) 

where  k  =  0.41  and  o)w  is  the  non-dimensional  wall  vorticity. 

In  the  numerical  computation  (T) ^  is  used  till  it  becomes  equal  to 
(T)q,  and  then  for  the  rest  of  the  shear  layer  the  outer  expression  has 
been  used. 


3.1.2  One-Equation  Modeling 

The  one-equation  model  is  basically  due  to  Glushko  [9  ]  in  the  form 
presented  by  Beckwith  and  Bushnell  [16]  and  recently  used  by  Coakley  and 
Bergmann  [14],  We  first  non-dimens ionalize  Eq.  (22)  as  follows: 


T  *  ,  E  *  -Ar  ,  P 

ij  u  2  ’  °  u  2  »  r 


-A-  P  ,  x  -  e  ,  D  -  -^r  5  (37) 

u  d  *  u  3  ’  u  3 


Thus  Eq.  (22)  in  non-dimensional  form  becomes 


3E  3E  3E  _  ,  _  ,  1  „2r 

H  +  u9j  +  v^-p  +  D-  x  +  iT''E' 

o 


(38) 


The  modeling  of  various  unknowns  proceeds  as  follows: 


T  =  |^  H(r) 
o 


(39) 


where 


(e)*^£  1*2 

- -  =  R  AE1, 

v  o 


R  ■  reference  value  of  R 
r 

D  »  JL 

R  R 

r 


(40a) 
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The  integral  length  A  is  defined  as 


i  <0.2 

-  0.2  .  A  >  0.2 


The  function  H(R)  is  defined  as 


(40b) 


H(R)  -  R  ,  R  <  0.75 

-  R  -(R  -  0.75)2,  0.75  <  R  <  1.25 

1  ,  R  >  1.25  J  (40c) 

In  the  modeling  of  the  diffusion  and  dissipation  terms,  a  function  G(R) 
will  be  needed,  which  is  defined  as 

G(R)  =  iiiMI 
H(R) 

=*  k  ,  R  <  0.75 


kR 


R  -  (R  -  0.75)2 


,  0.75  <  R  <  1.25 


kR  ,  1.25  <  R  < 


0.75 

k 


kR  -  (kR  -  0.75) 2  ,  <  FF 


1  ,  R  > 


1.25 


(40d) 


Having  defined  H  and  G,  and  using  the  following  abbreviated  notation 

x. 

yi  "  L  '*1  ~  7 2 

u. 


(yi  -  y2  ■  y) 
ui  "  F  (ui  “  u*  u2  "  v) 
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we  have, 

(i)  Reynolds  stress: 


(11)  Dissipation: 


(111)  Diffusion: 


3u.  3u,  » 

TiJ  -  -I(V  ^  +  3  E  Sil 


x  -  ff  +  kGT)  , 


c  =  constant  . 


D  -  (kGT 

3yk  a),k 


In  the  numerical  solution  of  the  energy  equation,  we  have  neglected  the 
spatial  variation  of  G.  The  values  of  the  constants  appearing  in  (39) 


and  (42)  are: 


k  =  0.41  ,  a  =  0.2  ,  c  *  3.93 


Substituting  the  preceding  expressions  in  Eq.  (38)  and  transforming  to  the 
curvilinear  coordinates  (5»n)»  we  have 

Et  ♦  <Vc  -  W  ■  5T  °>2e  +  hh  +  °lEr,>  +  TQ  -  x  (45) 

8  e 

where 


e  1  +  kGR  T 
o 


T1  "  T  +  kGRe(g22T£  “  812Tq) 

°!  "  0  +  kGRe(8llTn  ‘  g12TC) 

The  boundary  conditions  for  Eq.  (45)  are: 


n  -  nB  :  E  -  o 

n+n  :  E  +  o  ,  or  E  , 

00  oo 
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Eoo  *s  ProPortional  to  the  freestream  turbulence  level. 


3.1.3  Two-Equation  Modeling 

The  two-equation  model  used  in  the  present  research  is  due  to  Saffman 
[10]  and  Wilcox  [17].  Details  of  this  model  in  relation  to  the  present 
research  are  available  in  [ 18] .  In  non-dimensional  form,  any  equation  of 
the  present  model  can  be  represented  as 


xt+-m  %\  -  V„>  ‘  isr  <d2a  +  +  °iV 


where 


+  AA  +  B 


R 


(48) 


A  1  +  CXT 

T1  ‘  T  +  MRl(822TC  -  H2V 

°1  "  a  +  MRx(8llTn  "  gl2TcJ 

The  form  of  A,  B,  and  M  for  different  equations  are  as  follows: 
For  A  -  u  :  =  2  ,  A  =  2V2T,  M  =  4  , 

B  =  4(T  u  -  T  v  +  2  T  v) 
xx  y  yy  x  xy  y' 

For  A  -  E  :  =  1  ,  A  =  -(f1^2  ,  M  =  1  ,  B  =  2TQ 

For  A  =  <p  ;  =  1  ,  A  ■  -  Z^112  ,  M  =  1  ,  B  =  o 

In  this  model  the  non-dimensional  eddy  viscosity  is  given  by 

a  2E 


2* 


1/2 


aQ  -  0.3  ,  ^  =  5/3 


(49) 
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? 


and 


(Q  +  co2) 


(50) 


where  Q  has  already  been  defined  in  (34) . 

The  function  <p  which  has  its  own  rate  equation  is  the  turbulence 
vorticity-density  and  is  related  with  the  turbulence  dissipation  function 
X  as 


X  ~ 


Thus  in  the  turbulent  core  region,  the  behaviour  of  ij>  is  such  that 

,1/2 


4> 


1/2 


o  (a  E)J 
o  o 


kN 


where  N  is  the  normal  distance  from  the  wall  and  k  is  the  Karman  constant 
(-  0.41).  This  consideration  together  with  an  approximation  of  the  E 
equation  near  a  wall  yields  the  value  of  the  empirical  constant  aQ  which 
is  nearly  0.3.  Proceeding  in  the  same  manner,  an  approximation  of  the  in¬ 
equation  near  the  wall  region  yields  the  value  of  8^  which  is  nearly  5/3. 
The  boundary  conditions  for  $  are 


at  n 


2500  a  ,  .  , 

*  ■  C  R — >  <z-> 

o  r 


n  =  :  <t>  -*•  o 

where  Z r  is  a  roughness  length. 

3.1.4  Algebraic-Stress  Closure  Model 

It  has  been  well  known  for  quite  some  time  that  the  accuracy  of  physical 
representation  of  the  actual  turbulence  fields  depends  on  the  order  of 
closure  of  the  transport  equations  describing  turbulence.  Thus  the  model 
succeeding  the  two-equation  model  is  that  of  second-order  closure,  or  the 
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Reynolds-stress  closure.  This  type  of  higher-order  modeling  has  been  used 
quite  successfully  by  Launder,  et.  al.,  [19],  [20]  for  channel  and  boundary 
layer  type  flows.  Though  there  is  no  conceptual  difficulty  in  using  this 
model  along  with  the  mean  Navier-Stokes  equations,  the  discouragement  on 
its  use  stems  from  the  excessive  computational  requirements  in  the  form  of 
solving  three  or  six  additional  partial  differential  equations  for  two- 
or  three-dimensional  problems  respectively.  A  class  of  closure  which  can 
reasonably  be  thought  of  as  equivalent  to  the  second-order  closure  and 
called  the  "algebraic  stress  closure"  has  been  developed  by  Launder,  et. 
al.,  [21]  and  Warsi,  et.  al.,  [22]  and  is  described  below. 

The  basic  equation  for  the  present  consideration  is  now  Eq.  (24). 
Following  the  analysis  in  Warsi  [22],  the  Reynolds  stress  term  can 
be  approximated,  on  using  Eq.  (24),  as 


Tij  '  3  6  5ij  +  -  +  (Pij  3  P  6ij) 


(51) 


where  all  quantities  have  been  defined  earlier,  and  yq  =  0.6  ,  d^  =  0.8. 
The  accuracy  of  representation  (51)  depends  on  the  plausible  assumption 
that  the  substantive  rate  of  change  of  /e  in  following  the  mean  motion 
is  small  in  comparison  with  the  other  terms,  viz.. 


-  d  /.Til.  tii  de 

*  3F  <-3  «  5F 

e  e 


(52) 


Equations  (51)  form  a  system  of  non-linear  algebraic  equations  for 
the  determination  of  various  Reynolds  stresses  to  be  used  in  a  two- 
equation  model.  However,  assuming  the  Reynolds  stresses  in  P  to  be 
available  from  the  previous  iteration,  the  equations  can  be  solved  as 
linear  simultaneous  algebraic  equations  for  t^.  The  results  presented 
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K  ] 

I  i 

i  ' 

I ; 


below  are  for  the  case  of  two-dimensional  mean  flow  in  the  vorticity— stream 
function  formulation.  The  closure  for  the  unknown  terms  for  the  two- 
equation  model  is  due  to  Launder  [20].  All  quantities  appearing  below 
are  non-dimensional. 

Defining  the  quantities 


Su  •  *rt  ’  S12  -  *xy  •  S22  ‘  *„ 


3.1.9  Collection  of  Various  Formulae  and  Approximations 
All  the  formulae  listed  below  are  non-dimensional, 

(a)  Friction  Velocity: 


9  1  W 

u*2-ir 

o 


(58a) 


(b)  Local  Skin  Friction  Coefficient: 


!KI 

R  U  2 

o  p 


(58b) 


where 


U  -  MAX  (U  ) 
P  s 


s  f  g  n 


(c)  Reynolds  Shear  Stress: 


U1  U2 


-  T(*  -  i|>  ) 

yy  xx 


(58c) 


(d)  Pressure  Coefficient: 

The  wall  pressure  pw  non-dimensionalized  by  pu^2  is  given  by 


pw(5)  "  Pw(5reP  +  R  U  1/2  (g12“5 

o  sref  g 


8ll%)dC 


Now,  the  pressure  coefficient  is  defined  as 


p„(  >  -  p 


1  2 
2pu» 


Cp  “  2[pw«)  -  PJ  . 


(58d) 
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(e)  Transition  Factor: 


In  this  paper  the  Cebeci-Smith  transition  factor  y  [23]  was 
multiply  the  eddy  viscosity  to  account  for  the  transition  region, 
dimensional  variables,  y  is  given  by 


1  tr 


1  -  Exp[ 


-R  *66  U  _1,66  S  _1'34(S  -  S  ) 

o _ g _ tr _ try  j-s  ds, 

s.  U  J 


1200 


tr  p 


used  to 
In  non- 


(58e) 


where  Str  denotes  the  surface  distance  of  the  location  of  transition. 

(f)  Boundary  Conditions  for  E  and  x  Very  Near  to  the  Wall: 

In  the  immediate  vicinity  of  a  wall  the  low  Reynolds  number  effects 
on  turbulence  quantities  are  quite  significant.  Though  the  generated 
coordinate  mesh  is  very  fine  near  the  wall  region  for  all  cases  considered 
here,  the  viscous  resolution  still  demands  much  finer  meshes.  To  circumvent 
this  difficulty  we  have  used  approximate  analytic  expressions  for  E  and  x 
for  the  first  mesh  point  off  the  wall.  The  form  of  E  used  here  is  based 
on  the  analysis  of  Launder,  et.  al. ,  [21]  for  the  wall-viscous  effects. 

In  non-dimensional  variables  the  expression  for  E  is 


E  -  .0275 1 0)^1  2N2  ,  as  N  +  o  .  (58f) 

To  obtain  the  form  of  x  as  N  +  o,  we  use  Townsend's  expression  for  the 
dissipation  function  for  equilibrium  boundary  layers  [24],  which,  for  near 
wall  regions  is 

=  0.4E3/2 
X  N 

Using  (58f),  we  get 


X  =  .001824 1  u>w  |  3N2  ,  as  N  -*■  o  .  (58g) 

Note  that  in  principle  E  -*■  o  at  the  wall  in  such  a  way  that  x  00 .  However, 
as  observed  by  Daly  and  Harlow  [25],  the  tendency  of  x  to  become  infinite 
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4.  Numerical  Solution  of  the  Transport  Equation 

This  section  deals  with  the  computational  aspects  of  the  equations 
governing  the  mean  turbulent  flow  fields.  Besides  the  method  of  numerical 
solution,  we  shall  also  discuss  the  method  of  generation  of  initial  data 
for  the  solution  of  the  energy  and  the  dissipation  equations,  and  the  wall 
boundary  condition  for  the  vorticity  equation. 


4.1  Method  of  Solution: 

All  the  partial  differential  equations  are  first  discretized  by  using 
a  first  order  backward  difference  for  the  time  derivative  and  the  second 
order  central  differences  for  the  £  and  q  derivatives.  The  resulting 
system  of  difference  equations  is  then  solved  by  iteration  using  the  method 
of  point-successive  overrelaxation  (SOR) .  Using  the  model  equation  (12), 
the  SOR  difference  approximation  amounts  to  having  the  implicit  expression 


‘m"  -  »£>  +  m±j  (59) 

where  p  is  an  iteration  counter,  (i,j)  denote  the  spatial  position, R 

i » J 

is  the  residual  term  containing  the  previous  time  solution  and  the  previous 
iterative  solutions  at  the  neighboring  points,  and  W  is  the  acceleration 
parameter.  For  details  of  difference  approximation  of  the  two-equation 
model  refer  to  [18].  Convergence  of  the  iteration  is  obtained  by  satisfying 
the  inequality 


(p+D  _  x(p) 


lxvp 


(p+1) 


i*  j  i»  j 


<  e 


(60) 


for  all  points  of  the  field. 

In  the  solution  of  the  vorticity  equation,  viz.,  equation  (30),  it  was 
found  that  the  upwind  differencing  of  the  term  enhances  the  stability 

of  the  difference  approximation  [26].  Thus,  for  only  the  vorticity  equation 
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the  upwind  differencing  for  the  term  \p 

n  5 


has  been  used,  which  is 


ip 

r\  e 


( (Vt.j 
( (Vi,j 


'"l.j  ■  “l-l.J* 
*“1+1,1  '  “l.J* 


•  *n*° 

-  *„  *  0 


(61) 


The  method  of  solution  for  turbulence  problems  can  now  be  summarized 
as  follows.  Having  developed  the  numerical  coordinates  for  a  given  body, 
all  the  metric  coefficients  are  generated  numerically.  First  the  equations 
for  the  laminar  flow,  in  the  present  case  the  vorticity  and  stream  function 
equations,  are  solved  starting  from  time  t  »  o.  Usually  the  laminar  solution 
is  allowed  to  develop  to  a  quasi-steady  state  form  and  then  turbulence 
equations  are  introduced.  In  no  case,  except  in  [27],  the  turbulence  was 
introduced  at  t  =  o.  In  the  process  of  obtaining  the  laminar  solution, 
the  initial  data  for  the  E  and  x  equations  (described  in  Sect.  4.2  below) 
are  also  developed  simultaneously.  Having  available  the  velocity  field 
and  the  initial  data  for  E  and  X ,  it  becomes  an  easy  task  to  switch  on  to 
the  simultaneous  solution  of  all  the  equations.  The  choice  of  the  turbulence 
model,  location  of  the  beginning  of  transition,  and  the  amount  of  free- 
stream  turbulence  are  the  inputs  for  a  problem  which  can  all  be  handled 
easily  by  the  computer  program  used  in  the  present  research.  The  computer 
program  also  calculates  the  necessary  boundary  layer  parameters,  e.g., 
displacement  and  the  momemtum  thickness,  external  velocity,  etc. 


4.2  Initial  Data  for  the  E-and  x  Equations: 

The  approximate  generation  of  initial  data  depends  on  some  well  known 
forms  for  the  eddy  viscosity  and  the  dissipation  function  e.  From  (32a), 
we  have 

vT  =  i2(Q)1/2 


25 


? 


<v>2 


(62) 


where  olq  »  0.3.  Using  the  non-dimensionalization  Introduced  In  Sect.  3, 
we  obtain  from  Eqs.  (62) 


E  =  —  A2Q 
a 

o 

x  =  a2q^2  (63) 

As  noted  In  Sect.  4.1,  Eqs.  (63)  are  used  while  solving  the  laminar 
equations.  The  generated  values  of  E  and  x  are  then  used  as  the  initial 
data  for  the  solution  of  the  turbulence  model  equations. 


4.3  Wall-Vorticity : 

The  wall-vorticity  is  obtained  by  using  a  combination  of  the  Israeli 
method  [28]  and  a  modification  proposed  in  [18],  The  value  of  wall-vorticity 
is  assumed  to  be  correct  when  it  yields  a  vanishing  tangential  velocity  at 
the  wall.  Denoting  the  index  i  for  the  position  along  a  wall  and  p  the 
iteration  counter,  the  iterative  expression  for  the  wall-vorticity  is 
given  by 


(0 


(P+1) 

i 


. .  (p_1 ) 
“i 


v 


(P-1)  V 
(t)i 


(p) 

(t)i 


where  v^ 


is  the  tangential  velocity  component  at  the  wall. 


(64) 
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5.  Discussion  of  Results 


The  mathematical  models  of  various  orders  developed  in  Section  3  have 
been  used  to  compute  the  two-dimensional  incompressible  mean  turbulent 
flows  past  arbitrary  finite-bodies.  A  computer  program  based  on  the  finite 
difference  method  as  discussed  in  Section  4  has  been  used  to  obtain  the  flow 
fields  for  various  body  shapes.  The  difference  from  one  body  shape  to 
another  is  only  in  the  input  of  the  metric  data  and  the  coordinates  for 
each  body  which  are  obtained  by  separate  programs,  e.g.,  [12],  The  following 
three  shapes  have  been  extensively  studied  and  their  results  are  discussed 
below. 

(a)  Elliptical  cylinder  at  zero  Incidence. 

(b)  NACA  663-OI8  airfoil  at  zero  incidence. 

(c)  Circular  cylinder. 

In  all  the  above  cases,  the  free  stream  is  normal  to  the  axis  of  the  cylinders. 
5.1  Elliptical  Cylinder 

The  calculations  for  this  case  are  based  on  the  input  data  from 
Schubauer's  [28]  experimental  set  up.  The  elliptical  section  is  of  minor 
axis  10.11  cm,  major  axis  29.92  cm  and  is  exposed  to  a  uniform  stream  of 
Reynolds  number  Rq  of  1.59  x  105  based  on  the  minor  axis.  The  free  stream 
turbulence  level  E  is  0.1084  x  10" 3 . 

OO 

The  characteristic  length  L  for  this  problem  is  the  minor  axis  which 
has  been  used  to  non-dimensionalize  all  lengths.  A  highly  contracted 
coordinate  system  is  used  near  the  wall  region  to  resolve  the  viscous 
effects  in  the  boundary  layer.  The  coordinate  plot  is  shown  in  Figure  2. 

The  transition  model,  Eq.  (58e) ,  was  introduced  in  the  region 

1.5  <_  —  1.8,  where  s  is  the  arc  length  along  the  surface  measured 

s  s 

from  the  forward  stagnation  point.  From  —  -  1.8  to  —  =  3.3048  the  full 

L  Li 


turbulence  model  has  been  used.  The  laminar  solution  was  continued  from 
2L 

tt=  o  to  a  time  —  and  then  the  zero-equation  model  was  used  with  the 

woo 

3L 

transition  effects  included  to  a  time  — .  Based  on  this  solution,  the  one- 

wco 

equation  model  was  used  to  conclude  the  solution  process.  As  is  seen  from 
the  velocity  vector  plots  in  Figure  3,  the  separation  occurs  in  the  far 
downstream  region  very  near  to  the  trailing  edge.  It  is  because  of  this 
reason  that  we  have  compared  the  computed  pressure  coefficient  with  Schubauer's 
[  8]  curve  B,  and  the  comparison  seems  satisfactory,  (Fig.  4). 

In  Figure  5  the  coefficient  of  skin  friction  c^  is  compared  with  the 
Schubauer's  data  as  reproduced  by  Cebeci  and  Smith  [23].  The  skin  friction 
coefficient  in  the  turbulent  part  of  the  flow  is  computed  by  the  formula 
XJ^c^  =  0.5  Eg,  where  Eg  is  the  turbulence  energy  at  the  edge  of  the  sublayer. 
Figure  6  shows  the  turbulence  energy  contours  around  the  ellipse. 

5.2  NACA663-018  Airfoil 

The  results  presented  for  this  case  are  for  zero  incidence  at  a  chord 
Reynolds  number  Rq  of  1.2  x  106.  In  this  case  too,  there  is  a  strong 
contraction  of  coordinates  near  the  surface  so  as  to  have  at  least  20  n-lines 
in  the  boundary  layer.  Figure  7  shows  the  plot  of  the  coordinate  lines 
around  the  airfoil.  As  seen  from  Fig.  7,  the  coordinates  are  generally 
non-orthogonal . 

Experimental  data  on  NACA663~018  section  is  available  in  the  papers  by 
Bursnall  and  Loftin  [29]  and  Gault  [30].  Numerical  computation  on  this 
section  has  been  performed  by  Briley  and  McDonald  [26]  particularly  with 
the  idea  of  performing  the  Navier-Stokes  analysis  on  the  separation  bubbles. 

In  Ref.  [26],  a  laminar  boundary  layer  solution  is  patched  with  the  Navier- 


Stokes  separation  bubble  analysis  which  is  finally  patched  with  the  turbulent 


boundary  layer  solution  downstream  of  the  bubble  region. 
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In  the  present  computations  the  Navier-Stokes  solution  was  generated 

for  the  whole  surface  both  for  the  laminar  and  turbulence  cases  without  any 

indication  of  a  separation  bubble  at  Rq  =  1.2  *  106.  The  velocity  vector 

Plot  in  Figure  8  shows  no  indication  of  a  flow  reversal  or  recirculation 

region  at  the  expected  position.  The  transition  model  was  introduced  at 

the  expected  position  of  the  transitional  bubble  which  is  nearly  —  «  0.72, 

L 

where  L  is  the  chord  length* 

The  computed  pressure  coefficient  for  this  case  is  shown  in  Figure 
9.  Figures  10  and  11  show  the  turbulence  energy  contours  and  the  stream¬ 
lines  respectively.  Only  zero  and  one-equation  models  were  used  in  this 
study  with  the  free  stream  turbulence  E  *  .3375  x  10-5. 

CO 

5.3  Circular  Cylinder 

The  turbulent  flow  past  a  circular  cylinder  has  been  considered  by 

the  authors  in  [27].  Details  of  computational  analysis  are  available  in 

[18].  In  the  present  case  the  two-equations  model  as  developed  by  Saffman 

and  Wilcox  [17]  has  been  used.  Results  of  this  computation  are  shown  in 

Figures  12-16  for  a  Reynolds  number  R  *  1.4  x  io6  based  on  the  radius  of 

o 

the  cylinder. 


29 


6.  Conclusions 

Solutions  of  the  averaged  Navier-Stokes  equations  along  with  the 
equations  of  turbulence  quantities  have  been  obtained  for  various  body 
shapes  in  non-orthogonal  curvilinear  body-fitted  coordinate  systems.  In 
principle,  the  developed  computer  program  is  capable  of  predicting  the 
averaged  flow  and  turbulence  fields  around  any  fintie  two-dimensional  body, 
using  either  the  zero,  one,  or  two-equations  turbulent  closure  model.  In 
addition,  an  algebraic  stress  closure  model  has  also  been  developed  which 
can  be  used  with  any  two-equation  model  for  predicting  the  flow  fields 
having  large  separation  and  recirculating  regions. 


Appendix  A 


Let  x^  be  the  Cartesian  coordinates  and  5* 
coordinates.  Then  the  metric  coefficients  are 

Slj  3^  3^ 

ij  =  iii  iil 

3xk  3xk 


the  general  curvilinear 


(A-l) 


(A-2) 


g  =  det  (g^) 

2  3(xl’x2’x3> 

J  =  g  =  - 

3(C1,C2,C3) 

Using  the  summation  convention  on  repeated  indices, 
of  the  second  kind  are  given  by 


(A-3) 

(A-4) 

(A-5) 

the  Christoffel  symbols 


In  particular 
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Transformed  Plane 
(Natural  Coordinates) 


Figure  1.  Physical  and  Transformed  Planes. 


Figure  2.  Generated  Coordinate  Plot  for  the  Schubauer's  Ellipse. 
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Figure  4. 


Pressure  Coefficient,  C  .  •  Data  from  Ref.  [8] 
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Figure  14.  Coefficient  of  Pressure  (C  *)  Referenced  to  a  Rear 

Stagnation  Pressure  of  1.0.  is  p  -  p^,  ^  .  The  Total 
Drag  Coefficient  C  is  .293  where  C  =  D/p  r  . 


